clear all; clc; close

%load EcoliHyperShockEpi.mat
load CompileHyperShockData.mat
expIndex = 7;

[numCells,~]=size(Widths{expIndex}); % numCells is number of cells


%      R = 0.000555556; % The average unperturbed growth rate
%      deltaT = R.*(t_vec{i});
%      L0 = exp(deltaT);

[~,numPoints]= size(Time{expIndex});

locationOfShock = 21;%32;
     
     %numPoints = locationOfShock; % when plotting only before the jump. 
 
% figure
% hold on
% 
% xlabel('Time (s)')
% ylabel('L/L0')
% title(ShockMags(expIndex))

     
    %% calculating the average length as function of time
     avg_width = zeros(numPoints,1);
     N = zeros(numPoints,1);
     for j=1:numCells
         if (~isnan(Widths{expIndex}(j,locationOfShock))) && (~(Widths{expIndex}(j,locationOfShock)==0))
           
           normalized_width = Widths{expIndex}(j,:)/Widths{expIndex}(j,locationOfShock);
           
           %plot(Time{expIndex},Lengths{expIndex}(j, :)/Lengths{expIndex}(j, locationOfShock))
           %size(Lengths{expIndex}(j, :))
             
           normalized_width(isnan(normalized_width)) = 10;
           

        %plot(deltaT,normalized_length)
        for k=1:numPoints
            if normalized_width(k) ~= 10
           
         avg_width(k) = (normalized_width(k) +  avg_width(k));
         N(k) = N(k) + 1;
            end
        end
        end
       % plot(deltaT,normalized_length)
     end
     %avg_length
     avg_width = avg_width ./ N;
    % n
     %N
     
%      figure
%      hold on
%      
%      plot(Time{expIndex}, avg_length)
%      xlabel('Time (s)')
%      ylabel('L/L0')
%      title(ShockMags(expIndex))
%      
%      fig2pretty
    

     
        %% calculating the standard deviation corresponding to the mean above
     std_width = zeros(numPoints,1);
     N = zeros(numPoints,1);
     for j=1:numCells
         if (~isnan(Widths{expIndex}(j,locationOfShock))) && (~(Widths{expIndex}(j,locationOfShock)==0))
           
           normalized_width = Widths{expIndex}(j,:)/Widths{expIndex}(j,locationOfShock);
           
           normalized_width(isnan(normalized_width)) = 10;
           

        %plot(deltaT,normalized_length)
        for k=1:numPoints
            if normalized_width(k) ~= 10
           
         std_width(k) = std_width(k) +  (normalized_width(k) -  avg_width(k))^2;
         N(k) = N(k) + 1;
            end
        end
         end
       % plot(deltaT,normalized_length)
     end
    
     std_width = sqrt(std_width ./ N);
    % n
     %N
     
%      figure
%      hold on
%      
%      plot(Time{expIndex}, std_length)
%      xlabel('Time (s)')
%      ylabel('std')
%      title(ShockMags(expIndex))
%      
%      fig2pretty
     
     
     
     figure
     hold on
     
     plot(Time{expIndex}, avg_width, 'b-')
     
     plot(Time{expIndex}, avg_width - 2*std_width, 'r-')
     
     plot(Time{expIndex}, avg_width + 2*std_width, 'r-')
     
     xlabel('Time (s)')
     ylabel('W/W0')
     title(ShockMags(expIndex))
     
     hLeg = legend('example');
     set(hLeg,'visible','off')
     
     fig2pretty
    
     
     %% Export arrays
     
     expName = ['LBto' num2str(ShockMags(expIndex))];
     
     writematrix(Time{expIndex},['time_'  expName '.csv']) 
     writematrix(avg_width,['avgWidth_' expName '.csv']) 
     writematrix(std_width,['stdWidth_' expName '.csv']) 
     
